%% Figure for the draft, second part
% Land value: shrunken versus unshrunken
% 2017-6-5
clc; clear all; close all;
workpath00 = pwd;

addpath('toolbox_xt');
addpath('toolbox_plot');

%% Figure option
nyear = 6;

%% Load Data and preliminary stuff
% load('results_shrunken.mat');
load('results_shrunken_s3.mat'); %price estimation
load('results_estX_shrunken_s3.mat'); %covariates estimation
% three models (for both shrinkage/no-shrinkage)
% M1, M2, M3
% M3 is our preferred specification
% results generated from main_new201_Xestpred_unshrunken.m

% two models to select
temp_model1 = est_us_m5; %price
temp_modelX1 = estX_M3_unsh; %covariates

temp_model2 = est_s_m5; %price
temp_modelX2 = estX_M3; %covariates

% other info
ncity = size(est.ZZ2,1);

temp_y = est.YY;
temp_x = est.XX;

xtval = est.xtval;
xtflag = est.xtflag;


%% Land value computation - model 1 (unshrunken)
YY = est.YY;
XX = est.XX;


m_XX = mean(XX);
m_XX2 = mean([XX(:,14).^2, XX(:,14).^3]);

temp_est = temp_model1;

temp_sig2 = temp_est.model.MSE; %residual error variance

tab_area = [];
tab_ua   = [];
tab_land = [];

tab_lv = [];
for j=1:1:ncity;
    
    cityid = xtval(j);
    temp_sel  = tract.combofips==cityid;
    temp_grid = tract.min_dist(temp_sel);
    temp_lon = tract.lon(temp_sel);
    temp_lat = tract.lat(temp_sel);
    temp_d2c = tract.d2c(temp_sel);
    
    temp_area = tract.landsqmi_ua(temp_sel); %urban area
    temp_area_n = temp_area/sum(temp_area); % normalized area
    
    temp_area2 = tract.landsqmi(temp_sel); %urban area
    temp_area2_n = temp_area2/sum(temp_area2); % normalized area
    
    temp_n = size(temp_grid,1);
    
    % contructing distance dependent covariates for prediction
    % predX: we estimate and predict outside of this m-file.
    temp_predX1 = temp_modelX1.predX(temp_sel,1:14);
    temp_predX2 = temp_modelX1.predX(temp_sel,15:16);
    
    % land value at tract center in time t
    % time-invariant part of the value
    temp_val_nv = temp_predX1*temp_est.bet + temp_predX2*temp_est.bet2 ...
        + temp_est.c10*log(temp_d2c) + temp_sig2/2; % value does not depend on time
    
    temp_val = zeros(temp_n, nyear);
    for tind = 1:1:nyear
        temp_val_t = temp_est.sp.spfun_t(temp_grid,tind);
        temp_val_t = temp_val_t(:,j);
        temp_val_t = exp(temp_val_t + temp_val_nv);
        temp_val(:, tind) = temp_val_t;
    end
    
    % [1] Integrated using urban area, average over time
    temp_val_ua = sum(temp_val.*repmat(temp_area_n,1,nyear),1);
    temp_val1 = nanmean(temp_val_ua); %average over years
    
    % [2] Integrated using land area, average over time
    temp_val_land = sum(temp_val.*repmat(temp_area2_n,1,nyear),1);
    temp_val2 = nanmean(temp_val_land);
    
    % [3] Total using urban area, average over time
    temp_val_ua_tot = sum(temp_val.*repmat(temp_area*640,1,nyear),1); %acre correction
    temp_val3 = nanmean(temp_val_ua_tot); %average over years
    
    % [4] Total using land area, average over time
    temp_val_land_tot = sum(temp_val.*repmat(temp_area2*640,1,nyear),1); % acre correction
    temp_val4 = nanmean(temp_val_land_tot); %average over years
    
    % [5] 1, 5, 10, 20, 40 miles from the city hall
    temp_grid = [0,1/2, 5,10,20,40]';
    mean_temp_d2c = repmat(mean(log(temp_d2c)), length(temp_grid),1);
    
    pred_y = [];
    for vind = 1:1:16
        eval(['new_betX = temp_modelX1.bet',num2str(vind), '(temp_modelX1.xtval1 == cityid,:);']);
        new_X = [ones(length(temp_grid),1), log(1+temp_grid), mean_temp_d2c];
        new_y = new_X*new_betX';
        pred_y = [pred_y, new_y];
    end
    temp_predX1 = pred_y(:,1:14);
    temp_predX2 = pred_y(:,15:16);

    % time-invariant part of the value
    temp_val_nv = temp_predX1*temp_est.bet + temp_predX2*temp_est.bet2 ...
        + temp_est.c10*mean_temp_d2c + temp_sig2/2; % value does not depend on time
    
    temp_val = zeros(length(temp_grid), nyear);
    for tind = 1:1:nyear
        temp_val_t = temp_est.sp.spfun_t(temp_grid,tind);
        temp_val_t = temp_val_t(:,j);
        temp_val_t = exp(temp_val_t + temp_val_nv);
        temp_val(:, tind) = temp_val_t;
    end
    temp_val5 = nanmean(temp_val,2)'; %year average
    
    % [6] Integrated urban land value over time
    temp_val6 = temp_val_ua;
    
    % [0] Simple average
    temp_val0 = nanmean(exp(temp_y(xtflag == xtval(j),:)));
    
    % [00] Simple average (weigted)
    temp_a = exp(temp_y(xtflag == xtval(j),:));
    temp_b = exp(temp_x(xtflag == xtval(j),14)); %log(size)
    temp_val00 = sum( (temp_a.*temp_b)/sum(temp_b)); %weighted average
    
    % store
    temp_tab = [temp_val0,temp_val00, temp_val1, temp_val2, temp_val3, temp_val4, temp_val5, temp_val6];
    tab_lv = [tab_lv; temp_tab];
    
    tab_area = [tab_area; [sum(temp_area), sum(temp_area2)]]; %ua and land area
    tab_ua = [tab_ua; temp_val_ua];
    tab_land = [tab_land; temp_val_land];
    
end
tab_head = {'combofips', 'cmsa', 'pmsa', 'n', 'simple', 'simple (weighted)', 'LV (UA)', 'LV (Land)', 'LV-total (UA)', 'LV-total (Land)', ...
    'LV (center)', 'LV(1/2 mile)', 'LV (5 mile)', 'LV (10 mile)', 'LV (20 mile)', 'LV (40 mile)', ...
    'LV-2005', 'LV-2006', 'LV-2007', 'LV-2008', 'LV-2009', 'LV-2010', ...
    };
% combofips information
tab_name = {};
for j=1:1:ncity
    temp_sel = name.pmsafips==xtval(j);
    
    tab_name = [tab_name; [name.pmsafips(temp_sel,:), name.cmsaname(temp_sel,:), name.pmsaname(temp_sel,:), name.ntot(temp_sel,:)]];
    
end
big_tab = [tab_head; [tab_name, num2cell(tab_lv)]];

big_tab1 = big_tab;
tab_lv1  = tab_lv;

%% Land value computation - model 2
YY = est.YY;
XX = est.XX;

m_XX = mean(XX);

temp_est = temp_model2;

temp_sig2 = temp_est.model.MSE; %residual error variance

tab_area = [];
tab_ua   = [];
tab_land = [];

tab_lv = [];
for j=1:1:ncity;
    
% city specific predictors
        %m_XX = big_m_XX(j,:);
        %m_XX2 = big_m_XX2(j,:);
        
%         temp_XX  = XX(xtflag == xtval(j),:); %covariates in this city
%         temp_XX2 = XX2(xtflag == xtval(j),:);
%         temp_DD  = DD(xtflag == xtval(j),:);
        temp_ZZ  = ZZ(xtflag == xtval(j),:);
        
        % temporary variables

    cityid = xtval(j);
    temp_sel  = tract.combofips==cityid;
    temp_grid = tract.min_dist(temp_sel);
    temp_lon = tract.lon(temp_sel);
    temp_lat = tract.lat(temp_sel);
    temp_d2c = tract.d2c(temp_sel);
    
    temp_area = tract.landsqmi_ua(temp_sel); %urban area
    temp_area_n = temp_area/sum(temp_area); % normalized area
    
    temp_area2 = tract.landsqmi(temp_sel); %urban area
    temp_area2_n = temp_area2/sum(temp_area2); % normalized area
    
    temp_n = size(temp_grid,1);
    
    % contructing distance dependent covariates for prediction
    % predX: we estimate and predict outside of this m-file.
    temp_predX1 = temp_modelX2.predX(temp_sel,1:14);
    temp_predX2 = temp_modelX2.predX(temp_sel,15:16);
    
    % land value at tract center in time t
    % time-invariant part of the value
    temp_val_nv = temp_predX1*temp_est.bet + temp_predX2*temp_est.bet2 ...
        + temp_est.c10*log(temp_d2c) + temp_sig2/2; % value does not depend on time
    
    temp_val = zeros(temp_n, nyear);
    for tind = 1:1:nyear
        temp_val_t = temp_est.sp.spfun_t(temp_grid,tind);
        temp_val_t = temp_val_t(:,j);
        temp_val_t = exp(temp_val_t + temp_val_nv);
        temp_val(:, tind) = temp_val_t;
    end
    
    % [1] Integrated using urban area, average over time
    temp_val_ua = sum(temp_val.*repmat(temp_area_n,1,nyear),1);
    temp_val1 = nanmean(temp_val_ua); %average over years
    
    % [2] Integrated using land area, average over time
    temp_val_land = sum(temp_val.*repmat(temp_area2_n,1,nyear),1);
    temp_val2 = nanmean(temp_val_land);
    
    % [3] Total using urban area, average over time
    temp_val_ua_tot = sum(temp_val.*repmat(temp_area*640,1,nyear),1); %acre correction
    temp_val3 = nanmean(temp_val_ua_tot); %average over years
    
    % [4] Total using land area, average over time
    temp_val_land_tot = sum(temp_val.*repmat(temp_area2*640,1,nyear),1); % acre correction
    temp_val4 = nanmean(temp_val_land_tot); %average over years
    
    % [5] 1, 5, 10, 20, 40 miles from the city hall
    temp_grid = [0,1/2, 5,10,20,40]';
    mean_temp_d2c = repmat(mean(log(temp_d2c)), length(temp_grid),1);

    % constructing distance dependent covariates for prediction
    n_temp_pred = size(temp_grid,1);
    temp_pred_grid = log(temp_grid+1);
    temp_pred_ZZ = repmat(temp_ZZ(1,2), n_temp_pred,1);

    % only works for M3 model here
    new_x = [ones(n_temp_pred,1), temp_pred_ZZ, ...
            temp_pred_grid, temp_pred_ZZ.*temp_pred_grid, ...
            mean_temp_d2c, ...
            ];
    new_z = {[ones(n_temp_pred,1), temp_pred_grid]};
    new_g = repmat(cityid, n_temp_pred,1);

    pred_y = [];
    for vind = 1:1:16
        eval(['temp_m = temp_modelX2.model', num2str(vind), ';']);
        new_y = predict(temp_m, new_x, new_z, new_g);
        pred_y = [pred_y, new_y];
    end
    temp_predX1 = pred_y(:,1:14);
    temp_predX2 = pred_y(:,15:16);


    % time-invariant part of the value
    temp_val_nv = temp_predX1*temp_est.bet + temp_predX2*temp_est.bet2 ...
        + temp_est.c10*mean_temp_d2c + temp_sig2/2; % value does not depend on time
    
    temp_val = zeros(length(temp_grid), nyear);
    for tind = 1:1:nyear
        temp_val_t = temp_est.sp.spfun_t(temp_grid,tind);
        temp_val_t = temp_val_t(:,j);
        temp_val_t = exp(temp_val_t + temp_val_nv);
        temp_val(:, tind) = temp_val_t;
    end
    temp_val5 = nanmean(temp_val,2)'; %year average
    
    % [6] Integrated urban land value over time
    temp_val6 = temp_val_ua;
    
    % [0] Simple average
    temp_val0 = nanmean(exp(temp_y(xtflag == xtval(j),:)));
    
    % [00] Simple average (weigted)
    temp_a = exp(temp_y(xtflag == xtval(j),:));
    temp_b = exp(temp_x(xtflag == xtval(j),14)); %log(size)
    temp_val00 = sum( (temp_a.*temp_b)/sum(temp_b)); %weighted average
    
    % store
    temp_tab = [temp_val0,temp_val00, temp_val1, temp_val2, temp_val3, temp_val4, temp_val5, temp_val6];
    tab_lv = [tab_lv; temp_tab];
    
    tab_area = [tab_area; [sum(temp_area), sum(temp_area2)]]; %ua and land area
    tab_ua = [tab_ua; temp_val_ua];
    tab_land = [tab_land; temp_val_land];
    
end
tab_head = {'combofips', 'cmsa', 'pmsa', 'n', 'simple', 'simple (weighted)', 'LV (UA)', 'LV (Land)', 'LV-total (UA)', 'LV-total (Land)', ...
    'LV (center)', 'LV(1/2 mile)', 'LV (5 mile)', 'LV (10 mile)', 'LV (20 mile)', 'LV (40 mile)', ...
    'LV-2005', 'LV-2006', 'LV-2007', 'LV-2008', 'LV-2009', 'LV-2010', ...
    };
% combofips information
tab_name = {};
for j=1:1:ncity
    temp_sel = name.pmsafips==xtval(j);
    
    tab_name = [tab_name; [name.pmsafips(temp_sel,:), name.cmsaname(temp_sel,:), name.pmsaname(temp_sel,:), name.ntot(temp_sel,:)]];
    
end
big_tab = [tab_head; [tab_name, num2cell(tab_lv)]];

big_tab2 = big_tab;
tab_lv2  = tab_lv;



%% Figure 2-a: Shrunken versus Unshrunken, integrated
% tab_lv2
% 1'simple', 2'simple (weighted)', 3'LV (UA)', 4'LV (Land)', 5'LV-total (UA)', 6'LV-total (Land)', ...
%     7'LV (center)', 8'LV (1/2 mile)', 9'LV (5 mile)', 10'LV (10 mile)', 11'LV (20 mile)', 12'LV (40 mile)', ...
%     13'LV-2005', 14'LV-2006', 15'LV-2007', 16'LV-2008', 17'LV-2009', 18'LV-2010'
% sum(~isnan(tab_lv1(:,3)))
% sum(~isnan(tab_lv1(:,13)))
% sum(~isnan(tab_lv1(:,14)))
% sum(~isnan(tab_lv1(:,15)))
% sum(~isnan(tab_lv1(:,16)))
% sum(~isnan(tab_lv1(:,17)))
% sum(~isnan(tab_lv1(:,18)))

ind = 3;
n_tot = cell2mat(tab_name(:,4));
sel = ~isnan(tab_lv1(:,ind));

v1 = log(tab_lv1(sel,ind));
v2 = log(tab_lv2(sel,ind));
temp_ZZ2 = est.ZZ2(sel,:);

workpath00 = pwd;
figpath = [workpath00, filesep, 'figure_draft_rev', filesep, 'fig02'];
chk_dir(figpath);

fig1 = figure;
setmyfig(fig1, [2, 0.2, 8, 6.5])

% plot(log(temp_ZZ2(:,2)),[temp_model1.alp_jt_mat(:,1), temp_model2.alp_jt_mat(:,1)],'*')
scatter(temp_ZZ2(:,2), v1,100,'d','filled', ...
    'MarkerFaceAlpha',4/8,'MarkerFaceColor',rgb('grey'));
hold on
scatter(temp_ZZ2(:,2),v2,60,'filled', ...
    'MarkerFaceAlpha',8/8,'MarkerFaceColor',rgb('black'));

ly = nanmean(v2,2);
lx = temp_ZZ2;
temp_xlim = xlim;
x1 = temp_xlim(1);
x2 = temp_xlim(2);
xpred = linspace(x1, x2, 100)';
ypred = [ones(100,1),xpred]*( (lx'*lx)\(lx'*ly) );
bet = (lx'*lx)\(lx'*ly);
disp('se for fig_lv_average');
mbet = fitlm(lx(:,2),ly);
mbet.Coefficients
se1 = mbet.Coefficients{2,1};
se2 = mbet.Coefficients{2,2};

plot(xpred,ypred,'--', 'linewidth', 4, 'color', rgb('orangered'));
plot(xpred(1),xpred(2), '*', 'Markersize', 0.0001, 'linewidth', 0.0001, 'color', rgb('orangered'));
hold off
xlabel('Log Urban Area (square miles)');
ylabel('Log Average Urban Land Values ($ per acre)');

set(gca, 'linewidth', 2, 'fontsize', 15)

% l = legend('Unshrunken', 'Shrunken', ...
%     [num2str(bet(1), '%6.2f'), ' + ', num2str(bet(2),'%6.2f'), '* log(Area)']);
% set(l,'fontsize', 20, 'location', 'northeast');

l = legend('Unshrunken', 'Shrunken', ...
    [num2str(bet(1), '%6.2f'), '  +  ', num2str(bet(2),'%6.2f'), '* log(Area)'], ...
    ['(', num2str(se1, '%6.2f'), ')', '     ', '(', num2str(se2, '%6.2f'), ')'] ...
    );
set(l,'fontsize', 17, 'location', 'northeast');


ylim([8,24])

cd(figpath)
saveas(gca, ['fig_lv_average.png']);
cd(workpath00)


%% Figure 2-b: Shrunken versus Unshrunken, central
% tab_lv2
% 1'simple', 2'simple (weighted)', 3'LV (UA)', 4'LV (Land)', 5'LV-total (UA)', 6'LV-total (Land)', ...
%     7'LV (center)', 8'LV (1 mile)', 9'LV (5 mile)', 10'LV (10 mile)', 11'LV (20 mile)', 12'LV (40 mile)', ...
%     13'LV-2005', 14'LV-2006', 15'LV-2007', 16'LV-2008', 17'LV-2009', 18'LV-2010'
% sum(~isnan(tab_lv1(:,3)))
% sum(~isnan(tab_lv1(:,13)))
% sum(~isnan(tab_lv1(:,14)))
% sum(~isnan(tab_lv1(:,15)))
% sum(~isnan(tab_lv1(:,16)))
% sum(~isnan(tab_lv1(:,17)))
% sum(~isnan(tab_lv1(:,18)))

ind = 7;
n_tot = cell2mat(tab_name(:,4));

sel = ~isnan(tab_lv1(:,ind));

v1 = log(tab_lv1(sel,ind));
v2 = log(tab_lv2(sel,ind));
temp_ZZ2 = est.ZZ2(sel,:);

workpath00 = pwd;
figpath = [workpath00, filesep, 'figure_draft_rev', filesep, 'fig02'];
chk_dir(figpath);

fig2 = figure;
setmyfig(fig2, [2, 0.2, 8, 6.5])

% plot(log(temp_ZZ2(:,2)),[temp_model1.alp_jt_mat(:,1), temp_model2.alp_jt_mat(:,1)],'*')
scatter(temp_ZZ2(:,2), v1,100,'d','filled', ...
    'MarkerFaceAlpha',4/8,'MarkerFaceColor',rgb('grey'));
hold on
scatter(temp_ZZ2(:,2),v2,60,'filled', ...
    'MarkerFaceAlpha',8/8,'MarkerFaceColor',rgb('black'));

ly = nanmean(v2,2);
lx = temp_ZZ2;
temp_xlim = xlim;
x1 = temp_xlim(1);
x2 = temp_xlim(2);
xpred = linspace(x1, x2, 100)';
ypred = [ones(100,1),xpred]*( (lx'*lx)\(lx'*ly) );
bet = (lx'*lx)\(lx'*ly);
disp('se for fig_lv_halfmile');
mbet = fitlm(lx(:,2),ly);
mbet.Coefficients
se1 = mbet.Coefficients{2,1};
se2 = mbet.Coefficients{2,2};

plot(xpred,ypred,'--', 'linewidth', 4, 'color', rgb('orangered'));
plot(1,1, '*', 'Markersize', 0.0001, 'linewidth', 0.0001, 'color', rgb('orangered'));
hold off
xlabel('Log Urban Area (square miles)');
ylabel('Log City Center Land Values ($ per acre)');

set(gca, 'linewidth', 2, 'fontsize', 15)
% l = legend('Unshrunken', 'Shrunken', ...
%     [num2str(bet(1), '%6.2f'), ' + ', num2str(bet(2),'%6.2f'), '* log(Area)']);
% set(l,'fontsize', 20, 'location', 'northeast');

l = legend('Unshrunken', 'Shrunken', ...
    [' ', num2str(bet(1), '%6.2f'), '   +  ', num2str(bet(2),'%6.2f'), '* log(Area)'], ...
    ['(', num2str(se1, '%6.2f'), ')', '     ', '(', num2str(se2, '%6.2f'), ')'] ...
    );
set(l,'fontsize', 17, 'location', 'northeast');

ylim([4,25])

cd(figpath)
saveas(gca, ['fig_lv_halfmile.png']);
cd(workpath00)

%% Figure 2-c: Shrunken versus Unshrunken, Gradient
% tab_lv2
% 1'simple', 2'simple (weighted)', 3'LV (UA)', 4'LV (Land)', 5'LV-total (UA)', 6'LV-total (Land)', ...
%     7'LV (center)', 8'LV (1/2 mile)', 9'LV (5 mile)', 10'LV (10 mile)', 11'LV (20 mile)', 12'LV (40 mile)', ...
%     13'LV-2005', 14'LV-2006', 15'LV-2007', 16'LV-2008', 17'LV-2009', 18'LV-2010'
% sum(~isnan(tab_lv1(:,3)))
% sum(~isnan(tab_lv1(:,13)))
% sum(~isnan(tab_lv1(:,14)))
% sum(~isnan(tab_lv1(:,15)))
% sum(~isnan(tab_lv1(:,16)))
% sum(~isnan(tab_lv1(:,17)))
% sum(~isnan(tab_lv1(:,18)))

ind1 = 8;
ind2 = 10;
n_tot = cell2mat(tab_name(:,4));
sel = ~isnan(tab_lv1(:,ind1))&~isnan(tab_lv1(:,ind2));

v1 = (log(tab_lv1(sel,ind2))-log(tab_lv1(sel,ind1)))/9.5;
v2 = (log(tab_lv2(sel,ind2))-log(tab_lv2(sel,ind1)))/9.5;
temp_ZZ2 = est.ZZ2(sel,:);

workpath00 = pwd;
figpath = [workpath00, filesep, 'figure_draft_rev', filesep, 'fig02'];
chk_dir(figpath);

fig3 = figure;
setmyfig(fig3, [2, 0.2, 8, 6.5])

scatter(temp_ZZ2(:,2), v1,100,'d','filled', ...
    'MarkerFaceAlpha',4/8,'MarkerFaceColor',rgb('grey'));
hold on
scatter(temp_ZZ2(:,2),v2,60,'filled', ...
    'MarkerFaceAlpha',8/8,'MarkerFaceColor',rgb('black'));

ly = nanmean(v2,2);
lx = temp_ZZ2;
temp_xlim = xlim;
x1 = temp_xlim(1);
x2 = temp_xlim(2);
xpred = linspace(x1, x2, 100)';
ypred = [ones(100,1),xpred]*( (lx'*lx)\(lx'*ly) );
bet = (lx'*lx)\(lx'*ly);
disp('se for fig_lv_gradient');
mbet = fitlm(lx(:,2),ly);
mbet.Coefficients
se1 = mbet.Coefficients{1,2};
se2 = mbet.Coefficients{2,2};


plot(xpred,ypred,'--', 'linewidth', 4, 'color', rgb('orangered'));
plot(xpred(1),ypred(1), '*', 'Markersize', 0.0001, 'linewidth', 0.0001, 'color', rgb('orangered'));
hold off
ylabel('Log Land Value Gradient ($ per mile)');
xlabel('Log Urban Area (square miles)');
set(gca, 'linewidth', 2, 'fontsize', 15)
% l = legend('Unshrunken', 'Shrunken', ...
%     [num2str(bet(1), '%6.2f'), ' + ', num2str(bet(2),'%6.2f'), '* log(Area)']);
% set(l,'fontsize', 20, 'location', 'northeast');

l = legend('Unshrunken', 'Shrunken', ...
    [' ', num2str(bet(1), '%6.2f'), '   - ', num2str(abs(bet(2)),'%6.2f'), '* log(Area)'], ...
    ['(', num2str(se1, '%6.2f'), ')', '   ', '(', num2str(se2, '%6.3f'), ')'] ...
    );
set(l,'fontsize', 17, 'location', 'northeast');
ylim([-1,1])


cd(figpath)
saveas(gca, ['fig_lv_gradient.png']);
cd(workpath00)


%% Figure 2-d: Shrunken versus Unshrunken, Gradient
% tab_lv2
% 1'simple', 2'simple (weighted)', 3'LV (UA)', 4'LV (Land)', 5'LV-total (UA)', 6'LV-total (Land)', ...
%     7'LV (center)', 8'LV (1/2 mile)', 9'LV (5 mile)', 10'LV (10 mile)', 11'LV (20 mile)', 12'LV (40 mile)', ...
%     13'LV-2005', 14'LV-2006', 15'LV-2007', 16'LV-2008', 17'LV-2009', 18'LV-2010'
% sum(~isnan(tab_lv1(:,3)))
% sum(~isnan(tab_lv1(:,13)))
% sum(~isnan(tab_lv1(:,14)))
% sum(~isnan(tab_lv1(:,15)))
% sum(~isnan(tab_lv1(:,16)))
% sum(~isnan(tab_lv1(:,17)))
% sum(~isnan(tab_lv1(:,18)))

ind1 = 8;
ind2 = 10;
n_tot = cell2mat(tab_name(:,4));
sel = ~isnan(tab_lv1(:,ind1))&~isnan(tab_lv1(:,ind2));

% v1 = (log(tab_lv1(sel,ind2))-log(tab_lv1(sel,ind1)))/10;
% v2 = (log(tab_lv2(sel,ind2))-log(tab_lv2(sel,ind1)))/10;

v1 = tab_lv1(sel,ind1)./tab_lv1(sel,ind2);
v2 = tab_lv2(sel,ind1)./tab_lv2(sel,ind2);

% v1 = log(tab_lv1(sel,ind1)./tab_lv1(sel,ind2));
% v2 = log(tab_lv2(sel,ind1)./tab_lv2(sel,ind2));

temp_ZZ2 = est.ZZ2(sel,:);

workpath00 = pwd;
figpath = [workpath00, filesep, 'figure_draft_rev', filesep, 'fig02'];
chk_dir(figpath);

fig4 = figure;
setmyfig(fig4, [2, 0.2, 8, 6.5])

scatter(temp_ZZ2(:,2), v1,100,'d','filled', ...
    'MarkerFaceAlpha',4/8,'MarkerFaceColor',rgb('grey'));
hold on
scatter(temp_ZZ2(:,2),v2,60,'filled', ...
    'MarkerFaceAlpha',8/8,'MarkerFaceColor',rgb('black'));

ly = nanmean(v2,2);
lx = temp_ZZ2;
temp_xlim = xlim;
x1 = temp_xlim(1);
x2 = temp_xlim(2);
xpred = linspace(x1, x2, 100)';
ypred = [ones(100,1),xpred]*( (lx'*lx)\(lx'*ly) );
bet = (lx'*lx)\(lx'*ly);
disp('se for fig_lv_ratio');
mbet = fitlm(lx(:,2),ly);
mbet.Coefficients
se1 = mbet.Coefficients{2,1};
se2 = mbet.Coefficients{2,2};

plot(xpred,ypred,'--', 'linewidth', 4, 'color', rgb('orangered'));
plot(xpred(1),ypred(1), '*', 'Markersize', 0.0001, 'linewidth', 0.0001, 'color', rgb('orangered'));
hold off
xlabel('Log Urban Area (square miles)');
ylabel('Ratio of City Center to 10-miles Distant Land Values');
set(gca, 'linewidth', 2, 'fontsize', 15)
% l = legend('Unshrunken', 'Shrunken', ...
%     [num2str(bet(1), '%6.2f'), ' + ', num2str(bet(2),'%6.2f'), '* log(Area)']);
% set(l,'fontsize', 20, 'location', 'northeast');

l = legend('Unshrunken', 'Shrunken', ...
    [num2str(bet(1), '%6.2f'), '  +  ', num2str(bet(2),'%6.2f'), '* log(Area)'], ...
    ['(', num2str(se1, '%6.2f'), ')', '    ', '(', num2str(se2, '%6.2f'), ')'] ...
    );
set(l,'fontsize', 17, 'location', 'northeast');

ylim([-1,17])


cd(figpath)
saveas(gca, ['fig_lv_ratio.png']);
cd(workpath00)

